\documentclass[a4paper]{article} 

\usepackage{listings}
\lstset{tabsize=4, %
	frame=shadowbox, %把代码用带有阴影的框圈起来
	commentstyle=\color{red!50!green!50!blue!50},%浅灰色的注释
	rulesepcolor=\color{red!20!green!20!blue!20},%代码块边框为淡青色
	keywordstyle=\color{blue!90}\bfseries, %代码关键字的颜色为蓝色，粗体
	showstringspaces=false,%不显示代码字符串中间的空格标记
	stringstyle=\ttfamily, % 代码字符串的特殊格式
	keepspaces=true, %
	breakindent=22pt, %
	numbers=left,%左侧显示行号
	stepnumber=1,%
	numberstyle=\tiny, %行号字体用小号
	basicstyle=\footnotesize, %
	showspaces=false, %
	flexiblecolumns=true, %
	breaklines=true, %对过长的代码自动换行
	breakautoindent=true,%
	breakindent=4em, %
	%escapebegin=\begin{CJK*}{GBK}{hei},escapeend=\end{CJK*},
	aboveskip=1em, %代码块边框
	fontadjust,
	captionpos=t,
	framextopmargin=2pt,framexbottommargin=2pt,abovecaptionskip=-3pt,belowcaptionskip=3pt,
	xleftmargin=4em,xrightmargin=4em, % 设定listing左右的空白
	texcl=true,
	% 设定中文冲突，断行，列模式，数学环境输入，listing数字的样式
	extendedchars=false,columns=flexible,mathescape=true
	% numbersep=-1em
}
\usepackage{cite}
\usepackage{fancyhdr}
\usepackage{lipsum}
\usepackage{extarrows}
\usepackage{ctex}
\usepackage{amsfonts,amssymb} 
\usepackage{amsmath} 
\usepackage{cases}
\usepackage{geometry}
\usepackage{graphicx}
\usepackage{listings}
\usepackage{xcolor}
\usepackage{framed}
\usepackage{float}
\usepackage{ulem}
\usepackage{booktabs}
\usepackage{amsthm}
\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{width=10cm,compat=1.9}
\geometry{a4paper,left=3cm,right=3cm,top=3cm,bottom=3cm}

\begin{document}
	\begin{center}
		{\heiti {\huge 第五章编程作业}}
	\end{center}
	
	\begin{center}
		{\large 求数2101\ 叶陈昊\ 3210106359}
	\end{center}

	\section{ProblemA}

	代码见\texttt{ProblemA.m}.在\texttt{Matlab}命令行窗口输入\texttt{a}，
	得到拟合多项式的系数:
	$$
	\begin{aligned}
		a=
		\left[
		\begin{array}
			{c}
			a_{0} \\
			a_{1}\\
			a_{2}
		\end{array}
		\right]
		=
		\left[
		\begin{array}
			{c}
		2.175719932241656\\
		2.670413385241769\\
		-0.238443935926774
		\end{array}
		\right]
	\end{aligned}
	$$
	
	从而，通过正规方程组求解，在权重函数$\rho =1$的条件下，
	得到的二次拟合多项式为：
	$$ p(x) = 2.17572+2.67041x-0.238444x^2. $$
	
	用\texttt{Matlab}画出的散点和拟合多项式如图\ref{A}所示.
	
	\begin{figure}[H]
		\centering
		\includegraphics[width=0.9\textwidth]{imgs/ProblemA.png}
		\caption{问题A的拟合图像}
		\label{A}
	\end{figure}
	
	\section{ProblemB}

	代码见\texttt{ProblemB.m}.在\texttt{Matlab}命令行窗口输入
	\texttt{x\_}，得到拟合多项式的系数:
	$$
	\begin{aligned}
		x\_=
		\left[
		\begin{array}
			{c}
			a_{0} \\
			a_{1}\\
			a_{2}
		\end{array}
		\right]
		=
		\left[
		\begin{array}
			{c}
		2.175719932241668\\
		2.670413385241761\\
		-0.238443935926773
		\end{array}
		\right]
	\end{aligned}
	$$
	
	从而，通过QR分解超定方程组的系数矩阵求解，在权重函数$\rho =1$的条件下，
	得到的二次拟合多项式为：
	$$ p(x) = 2.17572+2.67041x-0.238444x^2. $$
	
	用\texttt{Matlab}画出的散点和拟合多项式如图\ref{B}所示.

	\begin{figure}[H]
		\centering
		\includegraphics[width=0.9\textwidth]{imgs/ProblemB.png}
		\caption{问题B的拟合图像}
		\label{B}
	\end{figure}

	在\texttt{Matlab}命令行窗口输入
	\texttt{cond\_G}，得到Gram矩阵$G$的条件数
	$$ \mathrm{cond}_2(G) = 1.89809\times 10^4. $$
	输入\texttt{cond\_R1}，得到矩阵$R_1$的条件数
	$$ \mathrm{cond}_2(R_1) = 1.37771\times 10^2. $$
	因此，验证得到$G$的条件数比$R_1$的大很多.
	
	\section{ProblemC}
	\subsection{$\alpha_j=\beta_j\ln 2+\gamma_j,\ 
		 \beta_j,\gamma_j$为有理数，$j=0,1,\cdots,n.$}\label{betagamma}
	 
	证明分为以下3步.
	\subsubsection{$r_j:=\int_{0}^{1}\frac{x^j}{1+x}\mathrm{d}x=(-1)^j\ln 2
	 	+\sum_{i=1}^{j}(-1)^{i+j}\frac{1}{i},\ j=0,1,\cdots, n.$}\label{rj}
	
	首先，$r_0=\int_{0}^{1}\frac{1}{1+x}\mathrm{d}x=\ln(1+x)|_{0}^{1}=\ln 2.$
	
	其次，根据提示，$$\sum_{i=1}^{j}(-1)^{i-1}x^{i-1}=\frac{1-(-1)^jx^j}{1+x}.$$
	将上式在$[0,1]$上积分，得：
	\begin{align}
		\int_{0}^{1}\sum_{i=1}^{j}(-1)^{i-1}x^{i-1}&=
		\int_{0}^{1}\frac{1-(-1)^jx^j}{1+x}\nonumber\\
		\sum_{i=1}^{j}(-1)^{i-1}\left.\frac{x^i}{i}\right| _0^1&=
		r_0-(-1)^jr_j\nonumber\\
		\sum_{i=1}^{j}(-1)^{i-1}\frac{1}{i}&=\ln2-(-1)^jr_j\nonumber
	\end{align}
	故
	$$ r_j = (-1)^j(\ln 2-\sum_{i=1}^{j}(-1)^{i-1}\frac{1}{i})= (-1)^j\ln 2
			+\sum_{i=1}^{j}(-1)^{i+j}\frac{1}{i},\ j=0,1,\cdots, n. $$
		
	\subsubsection{Hilbert矩阵可逆，且其逆矩阵中元素为有理数.}\label{H}
	我们先证明如下引理.\\
	\textbf{引理}：设$H$为$n$阶Hilbert矩阵
	（即$h_{ij}=\frac{1}{i+j-1}, 1\le i,j\le  n$），
		$ b=[0,0,\cdots,0,1]\in \mathbb{F}^n $，
		则线性方程组$Hx=b$的解$x=[x_1,x_2,\cdots,x_n]^T$存在唯一，且$x_i\in\mathbb{Q},\ i=1,2,\cdots,n$.\\
	\textbf{引理的证明}：构造函数
	$$ R(t):=\frac{x_1}{t+1}+\frac{x_2}{t+2}+\cdots+\frac{x_n}{t+n}. $$
	则$Hx=b\Leftrightarrow R(0)=R(1)=\cdots=R(n-2)=0,R(n-1)=1.$\\
	通分$R(t)$，得：
	\begin{align}
		R(t)=\frac{x_1(t+2)\cdots(t+n)+\cdots+x_n(t+1)\cdots(t+n-1)}
		{(t+1)(t+2)\cdots(t+n)}=:\frac{P(t)}{Q(t)}.\label{R}
	\end{align}
	从而
	\begin{align}		
		\mathrm{deg}(P(t))\le n-1,\ P(0)=P(1)=\cdots=P(n-2)=0;\label{P}\\
		\ P(n-1)=Q(n-1)=n(n+1)\cdots(2n-1).\label{P(n-1)}
	\end{align}
	由(\ref{P})可以得到$$ P(t)=c\cdot t(t-1)\cdots(t-n+2),$$
	其中$c$待定.结合(\ref{P(n-1)})得
	$$ c=\frac{n(n+1)\cdots(2n-1)}{(n-1)(n-2)\cdots2\cdot1}
	=\frac{(2n-1)!}{(n-1)!(n-1)!}=(2n-1)\binom{2n-2}{n-1}\in\mathbb{Q}. $$
	从而$P(t)\in \mathbb{Q}[t].$
	根据(\ref{R})中定义，
	$P(-k) = x_k\cdot\prod_{1\le j\le n,j\neq k}(-k+j),\ 1\le k\le n.$
	由于$P(-k)\in\mathbb{Q}, 
	\prod_{1\le j\le n,j\neq k}(-k+j)\in\mathbb{Q}\backslash\{0\}$，
	从而$x_1,\cdots.x_n$被唯一确定（即$x$存在唯一），
	且$x_k\in\mathbb{Q}, 1\le k\le n.$	$\hfill\square$
	
	根据以上引理的证明过程，我们不难得知，
	将向量$b=e_n$换成任意的$e_j,j=1,2,\cdots, n-1$
	(其中$\{e_1,\cdots, e_n\}$为$\mathbb{R}^n$的标准正交基)，
	引理的结论都是成立的.
	
	下面证明\ref{H}.设$HB=I_n,B=[B_1\ B_2\ \cdots \ B_n]
	\in \mathbb{R}^{n\times n}, I_n = [e_1\ e_2\ \cdots \ e_n].$
	则$HB_k=e_k,1\le k\le n.$由引理，$B_k$存在唯一，且其中元素为有理数.
	从而$B$存在唯一，且其中元素为有理数.
	进一步地，$B=H^{-1}\in \mathbb{Q}^{n\times n}.$
	
	\subsubsection{证明原命题}
	记\ref{rj}中$r_j = s_j\ln 2 + t_j,j=0,1,\cdots,n.$
	则$s_j,t_j\in \mathbb{Q},j=0,1,\cdots,n.$
	记$H_{n+1}$为$n+1$阶Hilbert矩阵，
	$$\alpha:=[\alpha_0,\alpha_1,\cdots,\alpha_n]^T\in\mathbb{R}^{n+1},$$
	$$ s:=[s_0,s_1,\cdots,s_n]^T, t:=[t_0, t_1,\cdots, t_n]^T
	\in\mathbb{Q}^{n+1}, $$
	$$ r:= [r_0,r_1,\cdots, r_n]^T=s\ln 2+t \in \mathbb{R}^{n+1} .$$
	所求解的正规方程组为:
	\begin{align}
		H_{n+1}\alpha = r.\label{Halpha=r}
	\end{align}
	由\ref{H}，$H_{n+1}$可逆，设其逆矩阵为$B\in\mathbb{Q}^{(n+1)\times(n+1)}.$
	则
	\begin{align}
		\alpha = Br
		=B(s\ln 2+t)
		=Bs\ln 2+Bt
		=:\beta\ln2+\gamma\nonumber
	\end{align}
	其中$\beta = Bs = [\beta_0, \beta_1,\cdots,\beta_n]^T\in \mathbb{Q}^{n+1},
	\gamma = Bt = [\gamma_0, \gamma_1,\cdots,\gamma_n]^T\in \mathbb{Q}^{n+1}.$\\
	故$\alpha_j=\beta_j\ln 2+\gamma_j,\ \beta_j,\gamma_j\in \mathbb{Q},
	j=0,1,\cdots,n.$
	$\hfill\square$
	
	\subsection{计算机求解$\beta_j,\gamma_j,j=0,1,\cdots,n.$}
	由\ref{betagamma},要想求解$\beta_j$以及$\gamma_j$，只需求解以下两个线性方程组：
	\begin{align}
		 H_{n+1}\beta&=s\label{Hbeta=s} \\
		 H_{n+1}\gamma &= t\label{Hgamma=t}
	\end{align}
	利用Gauss消去法，在\texttt{ProblemC.ipynb}中给出的精确解（用有理数表示）如下：
	\begin{figure}[H]
		\centering
		\includegraphics[width=0.9\textwidth]{imgs/ProblemC2beta.png}
		\caption{$\beta$的精确解}
		\label{beta}
	\end{figure}
	\begin{figure}[H]
		\centering
		\includegraphics[width=0.9\textwidth]{imgs/ProblemC2gamma.png}
		\caption{$\gamma$的精确解}
		\label{gamma}
	\end{figure}
	结合$\alpha = \beta\ln 2+\gamma$，得到$\alpha$的精确解如下：
	\begin{figure}[H]
		\centering
		\includegraphics[width=0.9\textwidth]{imgs/ProblemC2alpha.png}
		\caption{$\alpha$的精确解}
		\label{P2alpha}
	\end{figure}

	\subsection{精确表达式求解$\alpha_j,j=0,1,\cdots,n.$}
	取$\ln 2$为机器精度，计算图\ref{P2alpha}中的各项数值如下：
	\begin{figure}[H]
		\centering
		\includegraphics[width=0.9\textwidth]{imgs/ProblemC3alpha.png}
		\caption{$\alpha$的近似小数解}
		\label{P3alpha}
	\end{figure}
	观察上表的每一列，可以看出，
	当$n$越来越大时，$\alpha_j$会逐渐靠近并收敛于$(-1)^j$.
	
	\subsection{取$\ln 2$的近似来求解$\alpha_j,j=0,1,\cdots,n.$}
	取$\ln 2\approx 0.69315$，计算$r\approx 0.69315s+t$
	并直接求解线性方程组(\ref{Halpha=r})，得到的解$\alpha$如下：
	\begin{figure}[H]
		\centering
		\includegraphics[width=0.9\textwidth]{imgs/ProblemC4alpha.png}
		\caption{$\alpha$的直接求解}
		\label{P4alpha}
	\end{figure}
	明显看到，$|\alpha_j|$在$n$越来越大时会超出1很多.
	
	为此，计算$H_{n+1}(1\le n\le 6)$的条件数如下：（见\texttt{ProblemC4.m}）
	\begin{align}
		\mathrm{cond}(H_2) &= 19.281470067903967,\nonumber\\
		\mathrm{cond}(H_3) &= 5.240567775860627\times10^2,\nonumber\\
		\mathrm{cond}(H_4) &= 1.551373873892904\times10^4,\nonumber\\
		\mathrm{cond}(H_5) &= 4.766072502410004\times10^5,\nonumber\\
		\mathrm{cond}(H_6) &= 1.495105864145339\times10^7,\nonumber\\
		\mathrm{cond}(H_7) &= 4.753673565648391\times10^8.\nonumber
	\end{align}
	可以看到，随着$n$越来越大，Hilbert矩阵的条件数迅速增大，当取$\ln2\approx0.69315$
	时，会给线性方程组(\ref{Halpha=r})的右端项$r$一个微小的扰动，此时直接求解这个方程组，
	由于系数矩阵$H_{n+1}$的高条件数，势必会带来极大的误差.因此从图\ref{P4alpha}可以看出，
	$\alpha_j$在$n$越来越大时明显不收敛于$(-1)^j$了.

	\subsection{Tikhonov正则化直接求解$\alpha_j,j=0,1,\cdots,n.$}
	本题中，Tikhonov正则化方程为（为方便辨别，列向量用黑体向量表示）：
	\begin{align}
		\alpha'\boldsymbol{\alpha} + H^*H\boldsymbol{\alpha} = H^*\boldsymbol{r}.
	\end{align}
	其中$\alpha'=10^{-10}$为正则化参数，$H$为$n+1$阶Hilbert矩阵.令$A:=\alpha'I_{n+1}+H^*H\in \mathbb{R}^{(n+1)\times(n+1)},
	 b := H^*r\in\mathbb{R}^{(n+1)}.$则我们需要求解的方程为：
	 \begin{align}
	 	A\alpha=b.
	 \end{align}
 	直接求解它，得到$\alpha$的解如下：
 	\begin{figure}[H]
 		\centering
 		\includegraphics[width=0.9\textwidth]{imgs/ProblemC5alpha.png}
 		\caption{Tikhonov正则化后$\alpha$的直接求解}
 		\label{P5alpha}
 	\end{figure}
 	从中可以看到，对于$\alpha_0$和$\alpha_1$而言，已经开始分别恢复其收敛性.
 	\subsection{附言}
 	本来想全部用\texttt{Matlab}做的，后来发现要求解精确方程（用有理数表示），我就转用
 	\texttt{Python}了.不过之前也写了一套\texttt{Matlab}的代码，也附在文件夹下了.
%	$$
%	\begin{array}{c|c|c|c|c|c|c}  
%		n &1 & 2 &3 &4 &5 &6
%		 \\  
%		\hline  
%		\beta_0 &10 &75  &516  &3405   &21918    &138775 \\  
%		\beta_1 &-18&-408&-5700&-63480 &-618870  &-5526864 \\  
%		\beta_2 &   &390 &13620&273630 &4161360  &53241300 \\  
%		\beta_3 &   &    &-8820&-413280&-10780560&-207100321 \\  
%		\beta_4 &   &    &     &202230 &11865420 &379964971 \\  
%		\beta_5 &   &    &     &       &-4665276 &-328592881 \\  
%		\beta_6 &   &    &     &       &         &107975868 \\ 
%	\end{array}  
	
	
\end{document}